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ABSTRACT 



Advection-dominated accretion flows (ADAFs) are known to be convectively 



o 
o 
o 

CN ' unstable for low values of the viscosity parameter a. Two-dimensional axisymmetric 

^ hydrodynamic simulations of such flows reveal a radial density proflle which is 

. significantly flatter than the p oc r~'^l'^ expected for a canonical ADAF. The 

modified density profile is the result of inward transport of angular momentum 
by axisymmetric convective eddies. We present three-dimensional hydrodynamic 
. ^ simulations of convective ADAFs which are free of the assumption of axisymmetry. 

\^ ' We find that the results are qualitatively and quantitatively similar to those obtained 

■ with two-dimensional simulations. In particular, the convective eddies are nearly 

■ axisymmetric and transport angular momentum inward. 

O " 

o 
o 
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. 1. Introduction 

■ Analytical and numerical calculations have shown that advection-dominated accretion flows 
^ • (ADAFs) experience strong convection, especially for low values of the viscosity parameter 
^ a (Narayan & Yi 1994, 1995; Igumenshchev, Chen & Abramowicz 1996; Igumenshchev Sz 

Abramowicz 1999, 2000; Stone, Pringle &: Begelman 1999). The numerical work reveals 
that convection significantly modifies the structure of ADAFs (Stone et al. 1999; Narayan, 
Igumenshchev &: Abramowicz 2000, hereafter NIA). The change in the structure is understood to 
be a consequence of inward transport of angular momentum by convective eddies (NIA; Quataert 
& Gruzinov 2000, hereafter QG; Ryu & Goodman 1992; Stone & Balbus 1996). 
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All the numerical simulations done so far have corresponded to axisymmetric two-dimensional 
flows. NIA and QG emphasized that the inward transport of angular momentum seen in these 
simulations is a consequence of axisymmetry (cf. Stone & Balbus 1996) and that the results could 
in principle be qualitatively very different in three dimensions. We present in this paper the first 
fully three-dimensional numerical simulations of ADAFs. The details of our numerical technique 
are given in §2 and the results are described in §3. The paper concludes with a discussion in §4. 



2. Numerical Method 

We solve the non-relativistic time-dependent Navier-Stokes equations of an accretion flow in 
the gravitational potential of a central black hole: 

|+.V..- = 0, (2.1) 

= -VP + pV^ + Vn, (2.2) 

p^ = -PV-v + Q, (2.3) 

where p is the density, v is the velocity, P is the pressure, $ = —GM/r is the Newtonian potential 
of the central mass M, e is the specific internal energy, 11 is the viscous stress tensor (we include 
all components), and Q is the viscous dissipation function. Since we assume an ADAF, there 
is no radiative cooling. We adopt the ideal gas equation of state, P = (7 — l)pe, where 7 is 
the adiabatic index, and write the kinematic shear viscosity coefficient as = ac^/flK, where 
Cg = {P/pY^'^ is the isothermal sound speed and Q.k = {GM/r^)^^"^ is the Keplerian angular 
velocity. Equations (2.1)-(2.3) are independent of the black hole mass M when the radius r is 
scaled by the gravitational radius rg = 2GM/c^ and time is scaled by Vg/c. 

Each time step of the numerical computation is split into two sub-steps: hydrodynamical 
and viscous. The hydrodynamical sub-step is calculated by using the PPM algorithm (Colclla & 
Woodward 1984), and the viscous sub-step is solved by using a standard explicit algorithm for the 
viscous contribution to the equation of motion (2.2) and the energy equation (2.3). 

We solve equations (2.1)-(2.3) on a 3D Cartesian grid consisting of a sequence of six nested 
sub-grids. Each sub-grid A; is 26 x 26 x 52 in size and is uniform, extending over the range 
(0,/?^"*), (0,i?^"*), {-R^^R^*) in X, y, z, respectively. The innermost grid has i??"* = 26rg and 
each successive i?^"* is equal to 2i?^"*^, such that the outermost sub-grid extends to i?g"* = 832rg. 
Figure 1 illustrates the geometry of the grid. Due to numerical limitations we compute the flow 
only over a quarter of the full cubic volume and apply azimuthal periodic boundary conditions 
in the {x, z) and {y, z) planes. At the outer boundary R^^^ we use a free outflow boundary 
condition, and in the innermost sub-grid we assume that matter drains freely inside a sphere of 
radius r 2± 3rg. The inner boundary condition mimics the effect of the last stable orbit around a 
Schwarzschild black hole. 
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The time-step At for the numerical integration is chosen in accordance with the Courant 

condition for the hydrodynamical sub-step. For this At, the viscous integration is numerically 
stable because we choose sufficiently low values for the viscosity parameter a. We assume that 
mass is injected into the computational domain over an equatorial slender torus with radius 
rinj ~ SOOvg. The matter is injected at a fixed rate Minj and has rotational velocity (with respect 
to the z-axis) equal to 0.95 times the local Keplerian velocity. As a result of viscous interactions, 
part of the injected matter moves inwards and forms an accretion flow, while the rest leaves the 
computational domain through the outer boundary. Wc follow the time evolution of the flow on 
time scales longer than the viscous time scale at Vinj- This allows the flow to attain steady state 
in a time-averaged sense. Note that when p is scaled by Mi^j, the results are independent of the 
value of Minj- Therefore, the simulations depend on only two dimensionless parameters: a, 7. 

3. Results 

We have calculated four models (A, B, C, D) which differ in the choice of a and 7 (see 
Table 1). All the models show complicated time-dependent convective motions, analogous to those 
described in detail by Igunicnshchev ct al. (1996), Igumenshchev & Abramowicz (1999, 2000) 
and Stone et al. (1999). The amplitude of the variability is not large, and the average properties 
of the flow are well represented by time-averaging. To good accuracy, the time-averaged flow is 
axisymmctric in all the models. 

Figure 2 shows the time-averaged radial profiles in the equatorial plane (z = 0) of p, 0, and 
Cg for the four models. In the inner regions, r ^ lO^r^, the profiles are well approximated by 
power-laws: Q oc r~^/^, Cg oc r~^/^, p oc r~^, where the index f3 is ^ 0.5 for Models A and C 
and 0.8 for Models B and D. (For r > 2 x lO^r^, the flows are strongly affected by the outer 
boundary at r « 8 x lO^rg and are not power-law in character.) The scalings of J7 and Cg and that 
of p in Models A and C arc identical to those predicted by the analytical work of NIA and QG. 
The /? ~ 0.8 scaling of Models B and D is a little surprising, and suggests that (3 depends on the 
value of 7. A similar dependence of /3 on 7 was found also in the axisymmctric 2D simulations of 
Igumenshchev & Abramowicz (2000). 

Figure 3 shows a snapshot of the flow in the equatorial plane for Model C. The other three 
models are qualitatively similar. The most important feature is that the density and velocity 
fluctuations are almost axisymmetric, as seen by the fact that the eddies are significantly elongated 
in the azimuthal direction. This suggests that the present 3D results are not likely to be very 
different from those obtained with 2D axisymmetric simulations. The minimum length-scale of 
perturbations in the radial r and polar 9 directions depends on the viscosity parameter a, being 
larger for larger values of a. 

The left panel in Fig. 4 shows a snapshot meridional cross-section of Model C in the 
(y, 2:)-plane. The shading shows the distribution of the vorticity vector, u) = mtv. Regions with 
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a positive x-component of the vorticity uj^ are shown in grey, and regions with negative uJx are 
shown in white. It is seen that the flow is dominated by radial streams of inflowing and outflowing 
material which are sandwiched in the ^?-direction. The streams arc chaotic and move around so 
that they are washed out when a time average is taken. Models A and B are very similar to Model 
C, except that they have 2 — 3 inflowing streams at a given moment, whereas Model C has usually 
5 — 6 streams. 

Model D is qualitatively different from the other three models, as shown by the right panel 
in Fig. 4. What is plotted here is the time-averaged flow pattern and vorticity. We see that the 
model has quasi-stationary structures which survive even after averaging over about 100 periods 
of the Keplerian rotation at r = lOOr^ . The flow pattern is almost symmetrical with respect to the 
equatorial plane, and shows the presence of a stable and regular arrangement of circulation cells 
in the meridional plane. 

The direction of transport of angular momentum is given by the sign of the (r99)-component 
of the Reynolds stress tensor, {v'^v'^). NIA showed that this quantity is negative in 2D simulations 
of convective ADAFs, but noted that it was a consequence of the axisymmetry of the simulations 
and that real 3D convection might behave differently. In Figure 5 we show the radial dependence 
of the Reynolds stress averaged over 9 for the four 3D models described in this paper. We see 
that, except for a small boundary region, the stress is negative over the bulk of the flow in 
all the models. This confirms that the 2D results published in previous papers are physically 
meaningful. Figure 5 shows that the magnitude of the Reynolds stress is positively correlated with 
the values of both a and 7. Model D has a surprisingly weak stress. This may be connected to 
the quasi-stationary meridional circulation seen in this model. 



4. Discussion and Conclusions 

We have used numerical simulations to study three-dimensional hydrodynamic accretion flows 
with inefficient radiative cooling (ADAFs). The simulated flows have nearly Keplerian angular 
momentum on the outside, relatively weak viscosity, a = 0.01, 0.03, and an adiabatic index 7 
equal to either 5/3 or 4/3. We find that the convective motions are nearly axisymmetric (Fig. 
3). As a result, the basic properties of the flows are very similar to those found in previous 2D 
axisymmetric simulations (Igumenshchev et al. 1996; Igumenshchev & Abramowicz 1999, 2000; 
Stone et al. 1999; NIA). 

Our most important result is that convection transports angular momentum inward (Fig. 5). 
This property of convection was identifled by Ryu & Goodman (1992) in analytical work on thin 
accretion disks. Recently, NIA showed that convection in 2D axisymmetric simulations of ADAFs 
also moves angular momentum inward. They and QG argued, however, on the basis of previous 
work by Stone & Balbus (1996), that this result was a consequence of enforcing axisymmetry on 
the accreting gas. The present work shows that even in 3D, where the flow is not constrained to be 



-5- 



axisymmetric, convective motions do transport angular momentum inward. In fact, we find that 
2D and 3D simulations show good quantitative agreement in the estimate of the (r(/7)-component 
of the Reynolds stress tensor, and in the radial dependences of various fluid variables. 

In the models, the turbulent pressure due to convection, Pturb = (P''^'^)/3, where v' is the 
velocity fluctuation, is ^ 10~^ of the gas pressure P. Therefore, convection does not have a 
direct dynamical effect on the flow. However, as argued by NIA and QG, convection can have a 
very strong indirect effect on the structure of the flow if it causes inward transport of angular 
momentum. The accretion flow can then achieve a critically balanced state in which the convective 
transport almost completely cancels the outward transport of angular momentum by viscosity. 
This leads to a flat radial density proflle p cx r~^/^. Since the amount of transport possible with 
convection is limited, the viscosity parameter a must be less than a critical value acrit ^0-1 for 
the new structure to be realized (NIA). The value of a in an ADAF is not known at present. 
Model-fitting to the observed light curve and spectra of soft X-ray transients suggests that 
a ~ 0.2 - 0.3 (Esin, McClintock & Narayan 1997). Numerical MHD simulations in 3D should 
provide an independent estimate of a (assuming that viscosity in ADAFs is produced by magnetic 
stresses). If such computations confirm that a in ADAFs is large, then convection may turn out 
to have only limited dynamical effect on these flows. Conversely, if a < Ocrit) one may expect 
ADAFs to resemble the analytical solutions described by NIA and QG. 

Our results may be relevant also to the physics of rapidly rotating convective stars. In such 
stars convection may reduce or even suppress the outward transport of angular momentum by 
magnetic stresses (see Bisnovatyi-Kogan et al. 1979). 

An important property of convective accretion flows is the outward transport of energy 
by convection. We find that 2D and 3D simulations give similar results; the effective outward 
luminosity due to convection is ~ 10~^ — 10~^Mc^ for the values of a and 7 considered here, 
where M is the net mass accretion rate of the black hole. Very likely the energy will be radiated 
on the outside of the accretion flow. If so, the results have important implications for the spectra 
and luminosities of accreting black holes and neutron stars (Igumenshchcv &: Abramowicz 2000). 
Convection could also introduce variability on different time-scales (Igumenshchev &; Abramowicz 
1999). 

The simulated flows presented in this paper show some interesting dependence on the 
adiabatic index 7. Models with 7 = 4/3 have a steeper density profile (Fig. 2) than models with 
7 = 5/3. In addition. Model D (a = 0.01, 7 = 4/3) shows a quasi-stationary meridional circulation 
(Fig. 4), whereas other models have chaotic convective motions. The reason for these differences 
is unclear. 
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Table 1. Parameters of the models. 



Model a 7 



A 0.03 5/3 

B 0.03 4/3 

C 0.01 5/3 

D 0.01 4/3 
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Fig. 1. — A 2D illustration of the nested Cartesian grid used in the calculations. The example 
shown has three sub-grids, each of size 10 x 10. 
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Fig. 2. — Radial profiles of the equatorial density p, angular velocity ^/VLk and isothermal sound 
speed Cs/rVlK for Models A (solid lines), B (dotted lines), C (dashed lines) and D (long-dashed 
lines). The dot-dashed and dot-long-dashed lines in the top panel show power-law profiles, r~^'^ 
and r~^'^ , for comparison. 
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Fig. 3. — Snapshot view of the equatorial density (left) and velocity field (right) in Model C. The 

density contours are spaced with Alogp = 0.05. Arrows in the right panel show the dimensionless 
relative velocity (v ~ {v))/\{v)\, where {v) is the time- averaged velocity. The background rotation 
of the accretion flow is anti-clockwise. Therefore, material rotating clockwise/anti-clockwise in the 
figure has less/more specific angular momentum than the background. 
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Fig. 4. — Vorticity and momentum vectors in the meridional (y, z)-plane corresponding to a 
snapshot of Model C (left) and a time-average of Model D (right). Arrows correspond to rpv. 
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Fig. 5. — Radial profiles of the (r(^)-component of the Reynolds stress tensor {v'^v'^) for Models A 
(solid line), B (dotted line), C (dashed line) and D (long-dashed line). The stresses have been 
averaged over polar angle 9 and normalized to the square of the Keplerian velocity v\. 



